\section{Numerical Experiments}

\subsection{Numerical method for finding fixed points}

The results reported in this section depend on the calculation of positive
steady states for closed, single terminal-linkage class networks. We now
describe the algorithm used to solve for the associated fixed points.

The existence result of Theorem \ref{fp-exist-map} can 
be specialized to address closed systems. This requires
the redefinition of problem \eqref{convex-fix} and mapping \eqref{mapping}.

We redefine the optimization problem to be 
\begin{align}
	\underset{v\in\R^n}{\minimize} &\quad v^TD(\log(v)-\1)  \notag
\\                     \st &\quad YD v  = YA^Tr &:\ y \label{convex-fix-closed}
\\                         &\quad v \ge 0                     \notag
\end{align}.
We let $\Omega$ be the subset
Given an initial point $0< \hat v \in \Re^n$, and a small tolerance $\tau$, the algorithm proceeds as follows:
\begin{algorithmic}
  \STATE $r \gets \hat v$
  \STATE $v^\star \gets \hat v$<++>
  \WHILE{$\|YA_k v^\star\|_\infty>\tau$}
  \STATE $(v^\star,v^\star_0)\gets \text{unique solution of \eqref{convex-fix-closed}} $
	\STATE $(r,r_0) \gets \frac{1}{2}(r,r_0) +\frac{1}{2}(v^\star,v^\star_0)$
  \ENDWHILE
  \label{fixpoint-alg}.
\end{algorithmic}

The evaluation of the minimization step requires solving the linearly
constrained convex optimization problem \eqref{convex-fix}. Our implementation uses the PDCO
package \cite{pdco} for this purpose. 

Provided that the iterates satisfy $v^\star(\mu_k)>0$, and that
the minimization step is solved with sufficiently high accuracy, the optimality
conditions for \eqref{convex-fix} will imply that
\[\|Y^T\lambda-\log(v^\star(\mu_k))\|_\infty \leq \delta,\] for some small value of
$\delta$. Thus, if the iteration converges, at the fixed point \eqref{mak} will
be satisfied to a precision $\delta$ and \eqref{fb} to a precision $\tau$.

The aforementioned algorithm has been tested extensively on randomly generated
networks, with noteworthy success. Section \ref{scn:convergence} reports our
experiments.
\subsection{An example network}

\begin{figure} [h]
\centering
\includegraphics[width=3in]{./Graphics/paperNetwork.jpg}
\caption{Example network with two terminal-linkage classes}
\label{fig:network-small}
\end{figure}

\begin{figure} [htp]
\includegraphics[width=4.5in]{./Graphics/ConcentrationVsIterationExample.jpg}
\caption{Concentration convergence with iteration}
\label{fig:ConcentrationVsIteration}
\end{figure}

In this section we consider a toy network shown in Figure
\ref{fig:network-small}. For this network the number of complexes $n = 7$, the
number of terminal-linkage classes $l = 2$, and the stoichiometric subspace $S =
\operatorname{span}\{ A+E-C, C-A-D, B-C\}$ has dimension $s = 3$. Therefore,
the deficiency for this network is given by $\delta = 7 - 2 - 3 = 2$ and hence
neither of the deficiency 0-1 theorems can be applied to calculate equilibrium
points. However, intuition suggests that a non-zero steady state will exist
because of the weak reversibility of this network. We apply the fixed point formulation
described in Section 2 and use the numerical algorithm to solve for the
strictly positive steady state. Figure \ref{fig:ConcentrationVsIteration}
illustrates the convergence of the fixed point iterations to the steady state.\\
  
\begin{figure} [h]
  \includegraphics[width=4.5in]{./Graphics/EquilibriumVsTotalMassExample.jpg}
  \caption{Equilibrium dependence on total mass} \label{EquilibriumVsTotalMass}
\end{figure}

Figure \ref{EquilibriumVsTotalMass} illustrates the change in steady state as a
function of the \emph{total mass} in the system. The experiment shows that, as
the total mass is increased, species $A$, $B$ and $C$ adjust linearly to the
additional mass, while species $D$ and $E$ stay at the same levels.

This can be explained analytically by the fact that the vector $\mathbf{1}$ lies in range of $Y^{T}$. 

\subsection{Generating random mass-conserving chemical reaction networks with a
prescribed number of terminal-linkage classes}

We now describe the sampling scheme used to generate networks for the numerical
experiments in this section.  The scheme generates networks with
a prescribed number of strongly connected components, number of complexes, and number of species.

To generate a directed graph with $m$ nodes and $l$ strongly connected
components, we iteratively generate Erdos Reyni\footnote{An Erdos Reyni graph
is a directed unweighted graph. Each edge is included with probability $p$ and
all edges are sampled iid.}  graphs until we sample one with the required
number of strongly connected components. Once the graph has been sampled, the
edges are weighted by assigning independent and uniformly distributed weights
in the range $(0,10]$. 

To generate the stoichiometry we use a parameter $r$ that defines the maximum
number of species that will form each complex. For each complex $j$ we first
sample the number of species that will participate and then sample the subset.
All draws are done uniformly and independently.  Finally we assign positive
values independently to each participating specie, by using the absolute value
of a standard normal unit variance distribution. To ensure mass balance we
normalize the sum of the stoichiometry of the species that participate in a
complex to one. 

\subsection{Convergence of the fixed point algorithm on large single
terminal-linkage and multiple terminal-linkage class networks}\label{scn:convergence}

The fixed point iteration produces sequences that, up
to a small tolerance, satisfy \eqref{fb} at all iterates. Ideally it also
monotonically reduces the infeasibility with respect to \eqref{mak} until
convergence. Our extensive numerical experiments indicate that this is in fact
the behavior.

\begin{figure} [h] \centering \includegraphics[scale=0.3]{./Graphics/InfeasibilityVsIteration.jpg}
  \caption{Typical infeasibility of algorithm sequence, for network with a
  single terminal-linkage class} \label{fig:typical-infeas-single} \end{figure}

\begin{figure} [h]
  \centering
\includegraphics[scale=0.3]{./Graphics/InfeasibilityVsIterationMultiple.jpg}
\caption{Typical infeasibility of algorithm sequence, for network with two terminal-linkage classes}
\label{fig:typical-infeas-multiple}
\end{figure}

Figure \ref{fig:typical-infeas-single} displays the sequence of infeasibilities
$\|YA_kv_k\|_\infty$ generated by solving for a fixed point of a single terminal-linkage
network with $50$ species and $500$ complexes where at most $10$ species
participate in each complex. Figure \ref{fig:typical-infeas-multiple} displays
the analogous sequence for a network of equal size and two terminal-linkage classes.  We
have observed this (apparently linear) convergence rate consistently over all
generated networks regardless of the number of terminal-linkage classes.


\begin{figure} [h] \centering
  \includegraphics[scale=0.3]{./Graphics/SingleNetAvgIterationsVsNetSize.jpg} \caption{Average number
  of iterations for single terminal-linkage class networks}
  \label{fig:iteration-count-simple} \end{figure}

\begin{figure} [h] \centering
  \includegraphics[scale=0.3]{./Graphics/MultipleNetAvgIterationsVsNetSize.jpg} \caption{Average
  number of iterations for networks with two terminal-linkage classes}
  \label{fig:iteration-count-multiple} \end{figure}

We have also investigated the number of iterations necessary to
converge on networks of different sizes, with one and two terminal-linkage classes.
Figures \ref{fig:iteration-count-simple} and \ref{fig:iteration-count-multiple}
display the average number of iterations necessary for convergence for network
sizes ranging from $100$ to $5000$ complexes. Notably the average number of iterations
increases less than $10\%$ as the network size grows fifty-fold. (The average
was taken over $20$ instances per network size.)
This is typical of the interior-point optimization approach used by PDCO.

Our future work entails proving theoretical results on the existence of
positive equilibria for chemical reaction networks with
multiple terminal-linkage classes. However, our comprehensive numerical
experiments seem to indicate that even for networks with more than one
terminal-linkage class, there exists at least one positive fixed point.

Finally, while working on this manuscript, we were made aware of related work by Deng et. al. \cite{Deng}. Extending the work by Feinberg and Horn, they prove that weak reversibility is necessary and sufficient for existence of positive equilibria. While their proof is for closed systems with $b=0$, their extension using the $0$ complex is not clear as mentioned earlier for solving for any given admissible $b$. More importantly, our proof uses a less complicated convex optimization formulation which along with a fixed-point algorithm gives a method to calculate numerical solutions.\\

\hspace{-0.2in} \textbf{Acknowledgement}:  We gratefully acknowledge Anne Shiu and her student for helpful discussions and pointing out \cite{Deng}.


